/*  scfcpowc.cpp  by K.Tsuru  */
// function ID = 9106
/*****************************************************************************
SComplex class
It returns z^a = exp(a*log(z)).

Let z = x+i*y, a = c+i*d,
a*log(z) = (c+i*d)*log(x+i*y) 
= (c+i*d)*{(1/2)*log(x^2+y^2)+i*arctan(y/x)}
= {c*(1/2)*log(x^2+y^2)-d*arctan(y/x)}+i*{d*(1/2)*log(x^2+y^2)+c*arctan(y/x)}
= (c*r-d*t)+i*(d*r+c*t)
= p+i*q
where r = (1/2)*log(x^2+y^2), t = arctan(y/x) ), p = c*r-d*t, q = d*r+c*t.
Then
z^a = exp(p+i*q) = exp(p)*{cos(q)+i*sin(q)}
= polar(exp(p), q)

Notice the relation : arctan(y/x) = atan2(y, x) = arg(z)
*******************************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
#define UsesSimpleVersionInCpowC 1
SComplex Cpow(const SComplex& z, const SComplex& a){
    if( z.IsZero(9106) ) SNManager::SetError(SNManager::DOMAIN_ERR, "Cpow", 9106);// x == 0
    if(a.IsZero(9106)) return SComplex(1.0, 0.0);   // a == 0
#if UsesSimpleVersionInCpowC // 241(sec) faster a little
// simple version
    SComplex y, result;
    y = Clog(z);
    y *= a;
    result = Cexp(y);
    return result;
#else  // 244(sec)
    SDouble r = z.Norm();   // = x^2+y^2
    r = Log(r)/2.0;         // (1/2)*log(x^2+y^2)
    SDouble t = Arg(z), p, q;
    p = a.Real()*r - a.Imag()*t;
    p = Exp(p);
    q = a.Imag()*r + a.Real()*t;
    return Cpolar(p, q);
#endif
}